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Abstract 

We present numerical simulations in 3D settings where coronal rain phenom¬ 
ena take place in a magnetic configuration of a quadrupolar arcade system. 
Our simulation is a magnetohydrodynamic simulation including anisotropic 
thermal conduction, optically thin radiative losses, and parametrised heating 
as main thermodynamical features to construct a realistic arcade configura¬ 
tion from chromospheric to coronal heights. The plasma evaporation from 
chromospheric and transition region heights eventually causes localised run¬ 
away condensation events and we witness the formation of plasma blobs due 
to thermal instability, that evolve dynamically in the heated arcade part and 
move gradually downwards due to interchange type dynamics. Unlike earlier 
2.5D simulations, in this case there is no large scale prominence formation ob¬ 
served, but a continuous coronal rain develops which shows clear indications 
of Rayleigh-Taylor or interchange instability, that causes the denser plasma 
located above the transition region to fall down, as the system moves towards 
a more stable state. Linear stability analysis is used in the non-linear regime 
for gaining insight and giving a prediction of the system’s evolution. After 
the plasma blobs descend through interchange, they follow the magnetic field 
topology more closely in the lower coronal regions, where they are guided by 
the magnetic dips. 
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1. Introduction 


Coronal rain blobs and prominences consist of cool and dense plasma and 
can be formed due to thermal instability, as it was described in Parker (1953) 
and Field (1965). First, chromospheric and coronal matter gets heated due 
to a physical trigger (e.g. a solar flare) and gets redistributed to coronal 
heights, where it starts loosing energy and cools down through radiative 
losses and thermal conduction (Colgan et ah, 2008; Goldsmith, 1971; Levine 
& Withbroe, 1976; Townsend, 2009; Xia et ah, 2011). Thermal instabil¬ 
ity magnifies temperature perturbations and causes condensation in both 
hydro and magnetohydrodynamic regimes, isotropically and anisotropically, 
respectively (Hildner, 1974). Observational studies provided indications that 
prominences are composed of embedded cool plasma inside a flux rope with 
a sheared arcade on top (Chae et ah, 2001). Catastrophing cooling can take 
place within preformed flux ropes (Xia et ah, 2014b) or at the top of dynamic 
loop systems reducing the temperature of the plasma to transition region and 
chromospheric temperatures causing runaway events and the evacuation of 
the loop, as presented in Antolin et al. (2010); Mendoza et ah (2005); Miiller 
et al. (2004b); O’Shea et al. (2007); Schrijver (2001). 

These loops show strong emission in lines with characteristic temperatures 
T < lO^K (Muller et ah, 2003; Schrijver, 2001) and weak emission at lines 
corresponding to higher temperatures. The cool material is observable in 
EUV lines, and can either be a large scale prominence structure (Keppens & 
Xia, 2014; Xia et ah, 2014b) or form smaller blobs with sizes of the order of 
a few hundreds of kilometres (Antolin & Rouppe van der Voort, 2012; Fang 
et ah, 2013; Muller et ah, 2003; Schrijver, 2001; Tripathi et ah, 2009). These 
blobs appear bright in 304 A(EIT), while intensity variations appear to travel 
from the loop top towards the footpoints, suggesting cool temperatures and 
providing supporting evidence for heating-evaporation-condensation cycles 
(De Groof et ah, 2005, 2004; Mok et ah, 2008; Xia et ah, 2012). 

The importance of the magnetic field geometry in coronal rain and promi¬ 
nence events was highlighted in Kawaguchi (1970), as it guides the cool 
plasma falling from the corona back to the chromosphere (Leroy, 1972). 
Goronal rain can provide information about the magnetic fields support¬ 
ing it, as condensed blobs, with lifetimes of a few tens of minutes (Fang 
et ah, 2013), slide and deform along the magnetic field fines revealing their 
small-scale topology. They slide down with a wide range of speeds and with 
accelerations different from the one predicted by gravity alone, suggesting 
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the important influence of other (M)HD forces (Antolin & Verwichte, 2011). 
Both hot plasma upflows from chromospheric heights and multi-temperature 
downflows along magnetic field lines were observed to take place mainly in 
EUV wavebands (Beckers, 1962; Kamio et ah, 2011), with hot downflows 
only close to the loop top (Tripathi et ah, 2009). 

Magnetic field topologies used for the simulations of loop systems extend¬ 
ing to the solar corona have e.g. used force-free field approximations to ex¬ 
plore how heating affects EUV and soft X-ray views (Mok et ah, 2005). In a 
full flux rope configuration, Xia et al. (2014b) also demonstrated prominence 
formation and demonstrated its SDO/AIA appearance, in accord with ac¬ 
tual observations. Computational efforts have already successfully produced 
prominence condensations resulting from different heating functions (Antio- 
chos & Klimchuk, 1991; Xia et ah, 2012). Numerical studies (Antiochos et 
ah, 1999; Dahlburg et ah, 1998) have suggested that for the prominence to be 
formed through catastrophic cooling (Xia et ah, 2012) the heating function 
of the loop needs to be localised close to the footpoints of the loop at chromo¬ 
spheric heights, whereas the magnetic dip topology is not a necessity (Karpen 
et ah, 2001). Downflows of condensed matter were found to accelerate in the 
dipped part of the configuration in numerical simulations (Antiochos et ah, 
1999), while observational studies verify the accelerated sliding of the cool 
plasma towards both footpoints (Oliver et ah, 2014; Schrijver, 2001; Wang, 
1999) with propagation speeds in the range of 10 up to 200km/s (Kleint et ah, 
2014) and accelerations smaller than the effective gravitational one (Antolin 
& Rouppe van der Voort, 2012) in MHD setups. 

High resolution observations with instruments such as Hinode/Solar Opti¬ 
cal Telescope, CRISP/Swedish Solar Telescope, SDO and IRIS are available 
and provide insight in the dynamic fine-structures (Scullion et ah, 2014) of 
the solar corona. Numerical studies have also started to explore the physical 
conditions that play a key role in dynamic loop systems, such as the detailed 
magnetic topologies and the heating mechanisms, since they determine the 
temperatures, densities, and speeds of the condensed material (Eang et ah, 
2013), and influence relevant waves and instabilities driving the phenomena. 
With the already obtained knowledge about these physical conditions we 
want to further investigate the driving mechanism (s) for the creation (Xia et 
ah, 2011) and the evolution of the condensed blobs. 
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2. Setup: physical and computational aspects 

2.1. Governing equations and initial configuration 

We present numerical simulations in 3D settings where coronal rain phe¬ 
nomena occur in a quadrupolar magnetic environment. Simulations are done 
with the grid adaptive MPI-AMRVAC ^ code as described in Keppens et al. 
(2012) and updated in Forth et al. (2014). As initial conditions we take a 
2.5D, potential magnetic field, forming a quadrupolar arcade configuration 
with a dip, as demonstrated in figure 1, augmented with a ID stratified 
equilibrium, such that in ideal MHD we realise a force-balanced state. The 
velocity is set to zero throughout, and the magnetic field is taken as 
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where the r/—direction is the height in our simulation, Bq = 4:X 10“^T and the 
local angle between the (x, |/)-plane and the field lines is set to 7r/4 at x = 0, 
y = 2Lo/7r, determining BpQ and B^o uniquely, making Bpo — 

and Bzo = BpQ{e~^ — e~^). The density and pressure are calculated by 
the gravitational stratification. We initialise the ID arrays of density and 
pressure in the bottom boundary from hydrostatic equilibrium as follows: 
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The stratification gets a transition region temperature Tfr = 1.6 x lO^K 
at a chosen height htr = 0.27 x lO^m and adopts a constant vertical thermal 
conduction flux upwards, taken at k,(T)^ = 200Jm“^s“^. The temperature 
in the top boundary is later on restricted to be Ttop = 2 x 10®K. 

Nevertheless, we extend physically the ideal MHD equations by including 
non-ideal terms and study non-ideal effects including (a) optically thin ra¬ 
diative losses Q, (b) anisotropic thermal conduction and (c) a parametrised 


^https://gitorious.org/amrvac 
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heating function H, such that the evolution equation of the total energy 
becomes 

dE 

— + V • {Ev + ptotv - BB ■ v) = pg • V + V ■ {k ■ VT) - Q + H (5) 
ot 

where E = p/— 1) -\- pu^/2 + 5^/2 is the total energy density including 
internal, kinetic and magnetic energy densities accordingly for each term. 
We take 7 = 5/3 and set magnetic permeability pq = 1, making the total 
pressure given hy ptot = p+S^/2. The thermal conduction is field-aligned and 
it has the following form: K|| = The optically thin 

cooling uses a tabulated temperature dependence A(r) and as the equation 
Q oc njfA{T) suggests it is proportional to the squared hydrogen number 
density. We adopt the came cooling table A(r) as described and used in Xia 
et al. (2011), section 2.1 and figure 1 therein (solid line). 

The parametrized heating term has the following prescription 
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We have two different heating phases, as is evident by the heating function. 
First, we relax our system from begin time t = to to t = Eeiax (~ 51min) 
using only the background heating term H^g till a quasi-equilibrium state 
is reached. Afterwards, we add an extra localised heating term Hih, as 
visualised in the figure 1. For the background heating we have an ampli¬ 
tude oi Hq = 3 X 10“®Jm“^s“^, while for the localised heating we have 
Hi = 2 X 10 ^Jm ^s i.e. Hi ~ lO^i^o- The length and pressure units 
are taken as Lunit = lO^m and Punit = 0.03175Jm“^, while the scale height 
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Figure 1: Here we present in a 3D fashion the localised heating function and the magnetic 
field topology, to emphasise that the heating is applied at the footpoints of the enveloping 
arcade. Three different cross-sections are shown, namely the diagonal and the left and 
bottom boundary planes of our simulation domain. The image corresponds to a physical 
time of about 205min. 
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of the background heating is = bLunit and the localised heating is pre¬ 
scribed by a ramp function R{t) with a linear variation from zero to one, for 
a specific initial time and duration. For the rest of the heating parameters 
we take a 2 = and = 0.3, xi = -Xr = -A.2Lunit, Vh = 

\h = 0.25L^„jj. We assume fully ionised plasma with 10:1 H:He abundance. 

In figure 2 we present the temperature profile along a fine parallel to the 
y—axis passing through the middle of the simulation domain, i.e. (0,y,0), 
in three different times throughout our simulation. More specifically, we 
show a) t = to, b) t = Udax, c) t = tuob initiation and d) the static averaged 
temperature profile according to semi-empirical VAL-C model (Vernazza et 
ah, 1981), which only extends to low heights (< 2.6Mm) and is a plane- 
parallel model for quiet sun. We start at t = to from a temperature profile 
that includes chromosphere, transition region and corona and we apply a 
background heating with a scale length of 50 Mm till we reach a quasi¬ 
equilibrium state at t = treiax- The transition region is shifted to larger 
heights and heated. In physical units, the relaxation phase follows the 3D 
heated arcade for ss 51 min. Also, t = tuob initiation is actually ~ 205 min 
after t tj^dax • 

2.2. Boundary conditions and numerical aspects 

As boundary conditions we use ghost cells of 2 grid layers exterior to the 
physical domain everywhere in accord with the discretisation scheme. For 
numerical advancing of the partial differential equations we use a three-step 
Runge-Kutta scheme. As evaluation scheme to move from cell-centered to 
edge values we use a third-order-accurate limited reconstruction, as described 
in Cada & Torrilhon (2009). For fluxes we adopt an HLLCD scheme, which 
is a suitably mixed contact-resolving HLLC with fallback to TVDLF, hence 
D(iffused) hybrid scheme introduced in Meliani et al. (2008). As Courant 
parameter we use 0.9 of the CFL computed timestep. A diffusive approach 
is used for magnetic monopole control. 

For our simulation we typically impose boundary conditions on the prim¬ 
itive variables: (p, v^, Vy, v^, p, By, Bz). 

• Side boundaries, {xmin) ^max) ( and ^max) 

( ^d/y^Yiitt ^Lunit) 

For each physical quantity on the side surfaces in x—direction of the 
simulation box we choose the symmetric boundary conditions for (p, 
Vy, Vz, p, By, Bz) and asymmetric boundary conditions for (u^,, B^,), 
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Figure 2: Temperature profiles are shown at times a) t = to, b) t = treiax^ c) t = 
hiob initiations together with d) the static averaged temperature profile according to semi- 
empirical VAL-C model for comparison. 
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while taking periodic boundary conditions for the direction for all 
physical quantities. 

• Bottom boundary at: Umin = 0 

At the bottom, we use asymmetric boundary conditions for the three 
components of the velocity to set zero flow there. We flx the density 
and pressure as indicated by the gravitational stratification of the ini¬ 
tial state and the magnetic field is fixed as prescribed by the analytic 
relations given in Eq. (1), (2), (3). 

• Top boundary at: Umax = 

First, the pressure and density in the top layers is calculated giving a 
top temperature. This local temperature, but limited to Tfop is then 
used together with a discrete extrapolation of the stratification. The 
velocity is set asymmetric. For the magnetic field, we fix to its 
initial value, adopt a one-sided difference expression to set = 0. 
Finally, in order to set the normal component By we apply the central 
differenced condition V • 5 = 0. 

The simulation uses a base grid of 120 x 120 x 120 grid points, activating 
a 3-level mesh refinement based on mixed evaluation of weighted discrete 
second derivates involving density and the x— and y— magnetic field com¬ 
ponents p : Bx : By in a 0.6 : 0.2 : 0.2 ratio. The effective mesh size is 
480 X 480 X 480 corresponding to physical distances of 208km x 167km x 208km 
in a single highest resolution grid cell. 

2.3. Filter for analysis of thermodynamic evolution 

We use filters to analyse the presence of blobs in our simulations, which 
only start appearing after more than 205 minutes of physical time elapsed 
after extra heating is switched on. We define coronal rain blobs as cool 
and dense plasma suspended in the corona, after it condensed in-situ there. 
Condensation happens after a long period of gradual mass transport to the 
corona due to evaporation from transition region and chromospheric heights, 
as a consequence of the localised heating function that heats the footpoints 
of the large arcade, as shown in figure 1. 

Blobs are found as dense, cool matter well above transition region heights. 
Practically, the criteria used for the classification of a plasma element, ac¬ 
cording to our numerical resolution and grid size, as a blob element, are 
y > hTR{x,z,t) and p > Tpunit and T < O.IMK. More specifically, we use 
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criteria that take into account the number density, the temperature and the 
height of each grid cell to characterise a cell volume as containing blob mate¬ 
rial. First of all, we need to quantify the transition region height hmix, z, t). 
We do so in 3D by searching our simulation grid starting from the bottom 
looking for the height where the temperature and density gradient become 
maximum for the first time. If this technique fails, we set 9Mm as the tran¬ 
sition region height locally, which is typically an overestimate of the true 
transition region height found. After having quantified the transition region 
surface, we are able to calculate the total mass of the cool material and the 
number of the blobs, the latter by clustering of neighbouring plasma cells 
into larger parts. Furthermore, we define a spatial cutoff for blobs that are 
in agreement with the blob-filter criteria, but have a total volume smaller 
than eight individual grid cells, i.e. we only take into account blobs with a 
volume of at least eight grid cells in total. 

More specifically, the grid cells, that satisfy the blob-filter conditions and 
are above the set cutoff, are labelled as cool matter and kept in a list. With 
the cool material thereby mapped in our 3D grid, we can group the different 
grid cells into individual blobs at each saved snapshot and calculate their 
number. After this calculation is completed, we drop out groups with less 
than eight grid cells, remap the cool material and then recalculate the blob 
number. When all individual blobs are defined and mapped, we can collect 
their physical properties (such as the total mass of the blob, the average 
density) and calculate the centroid location of each blob with the temperature 
as weight. If the maximum density grid cell is close to the centroid and if their 
separation distance is larger than the resolution, we adopt as the centroid 
the grid cell with the maximum density value for that blob. 

In figures 3 and 4 we demonstrate the results of the analysis using the 
blob filter. Specifically, in figure 3 we show the total mass accumulated in the 
blobs (left panel) as well as the number of the blobs as the simulation evolves 
(right panel), while in figure 4 we present the distribution of the mass in blobs 
(left panel), as well as the distribution of the average temperature inside each 
blob (right panel) throughout the entire simulation. The minimum mass that 
can be contained in the smallest grid cell corresponding to the cutoff density 
criterion of the blob filter is about Arrimin = 1.1845 x lO^kg. This means that 
the peak in the blob mass distribution is about three orders of magnitude 
bigger than this minimum mass, i.e. nipeak oc 10^Arrimin, and is thus not 
influenced by resolution effects (also aided by the way we introduced a lower 
cut-off in blob size). In figure 3, we notice that there are two phases, namely, 
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Figure 3: The left plot shows the evolution of the mass contained in the blobs with time 
and the right one shows how the number of the blobs evolves with time. 
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Figure 4; The quantities shown on this figure are the temperature (right) and mass (left) 
distribution of all the blobs throughout the entire simulation. The blob material is rather 
cool with most of them having a temperature of 20,000K, while the range is from 10,000K 
to 100,000K. The vast majority of the blobs has a very small mass, close to the lower limit 
for their definition, while the distribution range is rather large. 


a build-up phase and a saturation phase. The build-up phase corresponds to 
mass accumulation at the middle of the domain, due to evaporation and blob 
formation. At some point, blobs obtain enough mass and are subject to the 
gravitational instability, so they start moving downwards. When they fall 
back into the transition region, our hlter stops taking them into account in 
the calculation of the total mass and the number of blobs. Thus, the episodic 
behaviour that we witness to take place in the left panel, after t = 230min and 
the saturation that we observe in the right panel of figure 3, after t = 235min, 
reflect the time variation of the blob creation and the blob loss, which sets 
up a mass recycling process. 
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3. Blob formation and evolution 


Due to evaporation from chromospheric and transition region heights, cool 
material starts concentrating at the top part of the large arcade preferentially 
within the heated loop part. As time passes, and more mass is accumulated 
at that specific region, at a height of about 4Mm, as shown in figure 5, the 
density starts rising and condensation centres make their appearance in our 
simulation. This formation of blobs containing cool material has been found 
and analysed in 2.5D simulations for a force-free arcade by Fang et al. (2013) 
and causes runaway events. In figure 5, the evolution of the mass transport 
and accumulation is presented in four different snapshots with a time interval 
of about 72min. The mass transport starts from the heated loop footpoints 
and moves upwards depositing mass, due to evaporation, on the loops above 
the magnetic dip of the arcade configuration, at a height of about 40Mm. 
The top left panel shows the moment when we activate the localised heating 
Hih {t = treiax) and reset the time to t = 0. The first blobs are formed at the 
time t ~ 205min and the last snapshot for our simulation corresponds to a 
time of about 259min. We have studied the phenomenon of the blobs for a 
time duration of about 54min, as we present further on. 

A continuous coronal rain develops, but here in a fully 3D fashion and 
the blob dynamics show clear indications of Rayleigh-Taylor or interchange 
instability. More specifically, as can be seen in figure 6, we observe verti¬ 
cal finger-like structures in the middle part of the velocity map (top panel), 
with opposite direction velocities appearing next to each other. Interesting 
features appear also close to the top of the simulation domain, with opposite 
velocities coexisting next to each other, indicative of a kind of stratifica¬ 
tion driven mixing occurring there. The magnetic field strength decreases 
exponentially with height, according to its analytic prescription that was 
presented in the previous section. In the second interesting region, at the 
top, the dominant component of the magnetic field is the third horizontal 
component, Bz- 

The temperature profile helps to locate the cooler blob as demonstrated 
in the bottom panel of figure 6, corresponding to the plane z — 0. Finger- 
like structures appear again at the same region as in the velocity profile. 
Interestingly, the first blobs make their appearance above the magnetic dip of 
our arcade configuration and this is the region from where on we start seeing 
the finger-like formations. From this perspective the condensed blobs seem 
to fall through the almost horizontal magnetic field lines, i.e the blobs look 
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Figure 5: Here we show the integrated view along the z— direction of the density at four 
different snapshots with a corresponding time interval of about 72min, to demonstrate the 
evolution and to show where the mass is accumulated due to the heating of the footpoints 
of the enveloping arcade. The time is annotated in minutes and the density is shown 
with the desaturated rainbow colorbar. In the top left panel, we show the t = Omin 
snapshot, when we activate the localised heating function Hih {t = Ueiax)- Later on, (top 
right and bottom left panel) the mass gets accumulated at heights around 40Mm. In the 
bottom right panel, the phenomenon of the blob runaway events has already started and 
condensations have made their appearance (at about 205 min blobs appear for the first 
time). 
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Figure 6: We present here the vertical component of the plasma velocity (top panel) and 
the temperature profile (bottom panel) at the time corresponding to 235min of physical 
time. Both Vy and T are shown in cross-section on the x — y plane corresponding to 2 ; = 0. 
We observe clear indications of Rayleigh-Taylor instability, seen as finger-like structures 
at heights ^ ^ 20 — 30Mm. 
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like they pass through the field lines, a scenario which would mean that the 
frozen-in condition gets violated. We will further analyse this phenomenon 
from different angles, to clarify the real blob motion later in this paper. 

4. SDO synthetic views 

We construct and present synthetic views from SDO/AIA in four different 
wavelengths 171 A, 193 A, 211 A, 304 A, that have their main emission con¬ 
tributions for plasma at temperatures of 0.8MK, 1.5 MK, 1.8MK, 0.08MK 
respectively. For the same time t = 235min as shown in figure 6 we present 
these four EUV views in figures 7 and 8, with the two first wavelengths 
(171 A, 193 A) corresponding to the first and second row of the figure 7 
and the two last ones (211 A, 304 A) corresponding to the first and second 
row of figure 8. We show two different sideway views corresponding to the 
z—integrated and the x—integrated view point for the left and right column 
of both figures accordingly, for one specific snapshot corresponding to 235min 
of physical time. The most informative view is the one that corresponds to 
the cool filter, i.e. the bottom row of figure 8. The blob dynamics and their 
structure becomes evident in this 304 A channel, allowing us to follow more 
closely the evolution of the position of the condensation centres with time. 
In animated views, we observe that mostly blobs form and move on the left 
part of the z—integrated view and on the right part of the x—integrated 
view. This means that (numerically) accumulated asymmetries have ren¬ 
dered one half of the heated arcade loops close to thermal instability. The 
cool material overall seems to follow the magnetic configuration during its 
evolution. In all filters, the footpoints of the enveloping heated large loop 
system appear enhanced in brightness, as the heating process is captured. At 
the 171A and 193A wavelengths (figure 7) the blob centres are enhanced and 
the phenomenon resembles coronal rain on these views (movie provided). 

5. Mass circulation 

As we will quantify further on, blobs form due to thermal instability, cre¬ 
ating condensation centres. After the blobs are formed, we will show further 
on that they are cospatial with regions liable to interchange instability. Its 
nonlinear evolution dominates the local dynamics and the blobs obtain an 
overall downward velocity. 
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Figure 7: SDO/AIA synoptic views in two different wavelengths showing the ^—integrated 
side view on the left and the x—integrated side view on the right of the whole simulation 
box. This shows a specific snapshot corresponding to the 235min of physical time in 171 
A (top panels) and 193 A (bottom panels) EUV channel. 
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Figure 8: SDO/AIA synoptic views in 211 A (top panels) and 304 A (bottom panels) wave¬ 
lengths showing the z—integrated side view on the left and the x—integrated side view on 
the right of the whole simulation box. This shows a specific snapshot corresponding to 
the 235min of physical time, and demonstrates that the cool blobs are mostly seen in the 
304AEUV channel. 
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time(min) 


Figure 9: The velocity component parallel to the local vector magnetic field is quantified 
here for five of the heaviest blobs in our experiment. We observe that there is a tendency 
for each blob in our sample to parallelise to the local magnetic field direction (parallel 
or antiparallel) later in their circulation through the corona as they approach the loop 
footpoints. 
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Figure 10: Here we present from two different lines of sight the magnetic field topology as 
well as the velocity vector field of the condensed material in four different snapshots. 
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More specifically, the blobs are very dynamic and as such they evolve 
changing their number, size and physical characteristics, splitting into smaller 
pieces or colliding forming larger cool and dense plasma elements. At the 
beginning, from the moment blobs start appearing in the simulation, they 
increase in mass and number, figure 3. The vast majority of the blobs have 
small masses, as shown in the left panel of figure 4, but overall the blobs 
have masses that spread out to a quite wide range of values varying from 
45 to 4500 10®g. As the time passes, they start splitting and small blobs 
start following more closely the magnetic field topology on the lower part 
of the configuration and specifically near the footpoints. Figure 9 provides 
evidence of this gradual alignment of the plasma clump motion with the 
local magnetic field direction, regardless their earlier directionality. At blob 
birth, blobs can have velocity components perpendicular to the magnetic field 
acquired when condensing (due to e.g. pressure difference). Note that we just 
quantify the centroid motion here. The blob growth during formation also 
influences motion of the blobs perpendicular to the magnetic field. Then they 
ultimately descend into the transition region with accelerating speeds. At 
time t ~ 225min, one part of the cool plasma follows the left footpoint of the 
large arcade to reach the transition region with increasing speed downwards 
while another part moves towards the centre of the configuration, just on top 
of the magnetic dip of the quadrupolar system, where a more complicated 
motion takes place. Specifically, overall we observe how blob features develop 
a rather circular motion on the top of the magnetic dip, that is accompanied 
by further blob splitting with some blobs falling along the field line topology 
towards both sides of the footpoints of the large arcade, while some parts 
have an apparent falling motion in the middle of the magnetic structure. This 
angular motion accelerates slightly. While this circular motion takes place 
new small blobs merge with the existing circulating plasma in that region. 
Another part of the cool plasma in the blobs that was found above the centre 
of the quadrupolar arcade system follows the topology on top of the small 
arcade and falls back to the transition region following their footpoints. The 
plasma that falls through the central part of the configuration, where the 
neighbouring footpoints with opposite polarities of the two small arcades 
come close, doesn’t show any noticeable acceleration. An impression of the 
blob dynamics is shown in Figure 10, which shows, from two perspectives, 
the field topology, and where each blob element has a fixed length vector 
representation with its direction defined by the flow, and colored by the flow 
magnitude. 
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The plasma-/?, representing the ratio of plasma pressure to magnetic pres¬ 
sure, in the lower part of the simulated region is smaller than unity, as indi¬ 
cated in figure 11. This is the region where the blobs get formed and evolve 
moving overall downwards. As the phenomenon progresses and the blobs 
move towards lower latitudes the plasma-/? decreases further in that region. 
More specifically the plasma-/? range of values is a) 8.4 x 10“^ < /? < 0.557 
for the blobs at all times b) 3.4 x 10“^ < /? < 1.7 for the physical domain 
above lOMm and c) 2.4 x 10“^ < /? < 121 for the whole physical domain. 
So we conclude that the blobs have low-beta values throughout the simula¬ 
tion. This observation is an indication that the magnetic field dominates. In 
order to accurately resolve and explain the plasma displacements and then 
link with the magnetic topology of our system, we need to develop particle 
tracing techniques on top of our MHD simulation to see how finid elements 
obey frozen-in conditions. In this paper we only present more indirect means 
to analyse the Lagrangian property and we leave this computationally de¬ 
manding particle tracing for future analysis. 



Figure 11: Plasma-/? profile on the slice that corresponds to the diagonal plane x = —z. 
The white curve indicates the region where /? = 1. Yellow represents high P values and blue 
low P values. We depict (left panel) one of the first snapshots that show the blob formation 
(?^214min), as well as, (right panel) the final snapshot of our simulation (?^259min), when 
there are still blobs circulating inside the coronal volume. Blobs are represented with 
white colour and they are projected on the 2D view. The plane x = —z has opacity 0.7 
so that it allows for blobs that are behind it to be seen in dark blue shades. 
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6. MHD stability analysis 

In the remainder of this paper, we will now make direct contact with 
linear MHD theory governing especially the instability criteria that may drive 
further nonlinear evolution. In general, the linearised MHD equations (about 
time independent physical quantities in equilibrium states p(f), .B(r), p(f)) 
accept solutions in normal mode fashion, according to Goedbloed & Poedts 
(2004): 

f(f, (11) 

where ^ is the plasma displacement vector. In ideal MHD, this leads to the 
eigenvalue problem: 

F(i) = ( 12 ) 

where F is the ideal MHD force operator. Both discrete and continuous 
eigenvalues are accepted and their set of {oJ^} form collectively the spectrum. 
For true ideal MHD conditions, the operator p~^F is self-adjoint so that the 
eigenvalues are real, which means the eigenvalues u) can be either real or 
imaginary, which respectively lead to waves and instabilities. 

For instabilities we have an exponential growth of the initial perturbation. 
However, our current setup is characterised by inclusion of non-ideal effects 
like thermal conduction and radiative losses. This makes that exponentially 
growing entropy modes can exist, or that magnetosonic wave modes can 
become overstable. Moreover, bulk flow also complicates the actual linear 
spectrum quantification (Goedbloed et ah, 2010). 

Still, we will now use earlier linear instability criteria for idealised setups, 
to quantify and study the first stages of the evolution of each perturbation. 
This will demonstrate that our simulation has thermal and interchange in¬ 
stabilities happening simultaneously. 

6.1. Instabilities and frequency eutoff 

The physical mechanism that creates the blobs and the complex dynam¬ 
ics that drives them around in the lower corona is related to two kinds of 
instabilities. When unstable, the system has a natural tendency to move 
towards more stable states through the mixing resulting from the developing 
instabilities. 

Since we have prominent radiative losses, as well as a strong magnetic held 
in the presence of gravity, we will consider: 

1. the thermal instability and 
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2. the interchange instability. 


The first one likely creates the condensation centres and thus the blobs and 
the second one can drive their evolution and circulation in the corona. More 
specifically, due to localised cooling and mass accumulation, the criteria for 
thermal instability can be fulfilled for condensation of cool material to take 
place in situ. This creates the observed blobs at first. After they are formed, 
they possess a higher density compared to the layers underneath them, so 
this configuration can become gravitationally unstable and the system needs 
to adjust to a new more stable state. In search of this new energetically 
economic configuration, the blobs start moving even in the region of strong 
magnetic field at the lower part of the simulation domain. The characteristic 
of the interchange instability is to do this with none to minimal field line 
bending k ■ B, where k denotes the wave vector when a plane wave type 
description is used. 

The spatial resolution of our simulation limits the wavelengths detectable 
in our data to (a multiple of) the grid cell size, Ax. Furthermore, data storage 
limits introduced a time cadence Atg ~ 85s, between successively stored full 
3D data, so the instabilities we can detect most easily have growth rates 
N in that grow sizeably in about 85s. Specifically for the Lagrangian 
displacement: 

e = , (13) 


we observe noticeable dynamic changes evolving from one saved frame to the 
next one when their displacement gets larger than our spatial resolution or 
when: 

Ax 


^NAts 


6 ■ 


(14) 


Initially, since blobs are created in situ, their Lagrangian displacement is 
small, so: 

|>1. (15) 

So, finally we get: 


eNAts 0.008Hz . (16) 

L\Lg 

We use this criterion further on to cut off lower growth rates artificially 
as they have not enough time to be determined by the system’s evolution 
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between two snapshots. Specifically, we are only interested in frequencies 
higher than 0.008Hz, which means that we are looking for negative squared 
eigenvalues in the region: < —6.4 x 10“^Hz^. These eigenvalues are the 

most important ones for determining the blobs’ formation and dynamical 
evolution at the temporal and spatial resolution that we have in this simula¬ 
tion. 


6.2. Thermal Instability: isochoric and isobaric criteria 

The prescribed localised heating function that heats up the footpoints 
of the external magnetic loop causes evaporation that gradually adds cool 
plasma from the transition region into the large loop of our main configu¬ 
ration. So the evaporated cool plasma increases the density in that region 
from about p = {Punit = 2.342 x 10“^^kg/m^) at t = 0 to p = punit 

at t ~ 205min. The density increase enters gradually the radiative loss and 
condensation centres are formed and blobs are created. 

Two thermal instability criteria were introduced, the isochoric by Parker 
(1953) and the isobaric one by Field (1965), that control thermodynamic 
evolution in astrophysical systems. From these criteria, we should be able 
to predict the condensation regions, i.e. regions where catastrophic cooling 
takes place. Although the criteria are only valid for gaseous uniform media, 
we will now use them to locate where blobs will be formed due to thermal 
instability 

According to Parker (1953) and also used for prominence onset in Xia et 
al. (2011) the isochoric criterion is: 


K 1 ar ar' 


(17) 


where k is the wavenumber of the perturbation that undergoes thermal 
instability, k, is the heat conduction coefficient, and R = nnUeRiT) = 
Qine/nn) is the radiative loss. 

Similarly, (see equation (8) in Xia et al. (2011)) the criterion for an isobaric 
thermal instability is given by: 

Cisobaric = P ~ ^ < 0, (18) 

where C = {nHneA{T) — H)/p = {R — H)/p = {Q{ne/nH) — H)/p is the 
generalised heat-loss function. We use the parallel conduction coefficient 
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value to quantify k (isotropic in a gaseous medium as analysed by Parker 
(1953)). 

As perturbation lengths (i.e. to set the wavenumber k) we use two different 
length scales that play an important role in the dynamic processes. 

1. On one hand, we adopt as perturbation wavelength the smallest de¬ 
tectable blob size A = 208 km, which corresponds actually to our spatial 
resolution. 

2. On the other hand, we can estimate the whole heated region size to 
be about A = 40Mm. This is done as follows. The temperature dis¬ 
tribution indicates that most blobs have an average temperature of 
about 20,000K, but in general blobs with temperatures from 10,000K 
to 100,000K are observed, with the vast majority of them concentrated 
in the range between 10,000K and 60,000K, as demonstrated in the 
second panel of figure 4. In order to estimate the size of the thermally 
unstable region, we use an isotemperature contour that corresponds to 
40,000K in our physical domain at the moment that we start observing 
for the first time the condensation phenomenon (205min). We esti¬ 
mate the size of the cool material region from the volume enveloping 
this isotemperature contour. The isosurface of 20,000K is present on 
the snapshot just before the first blob is formed, but vanishes on the 
snapshot that first captures the blob formation. This is in agreement 
with previous studies and is explained by the fact that cooling takes 
place in order for the condensations to happen, up to a point where 
cooling stops and inflows start converging on it, which results to a pres¬ 
sure increase that results to a localised increase in the temperature. In 
the same region the density is increasing strongly, forming condensation 
centres that gradually evolve into blobs. 

Using the isochoric thermal instability criterion, we conclude that the 
thermally unstable region is very well identified, as we can see on figure 12. 
Indeed, for the wavelength A = 40Mm, the condensation centres are captured 
and followed throughout the evolution of the phenomenon, as new high den¬ 
sity blobs keep emerging till the end of our simulation. 

The isochoric criterion is rather independent of the length scale of the 
perturbation in our analysis, as after using the two wavelengths that were 
estimated above, that differ about two orders of magnitude from each other, 
we do not observe large differences in the estimated unstable region with the 
isosurfaces’ method. This fact is in agreement with Xia et al. (2011), who 
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first showed the agreement with Parker’s results in a ID setup. So the value 
of the perturbation wavelength A and thus the perturbation wavenumber 
k = 2 tt/X, doesn’t affect the results of the thermally unstable analysis in 
highlighting the regions where the blob formation will take place. 

For the isobaric case though, we noticed a clear dependence on the length 
scale of the perturbation. Specifically, for the wavelength of A = 208km, the 
isobaric criterion suggests fully stable states throughout the simulation, as it 
gives Cisobaric > 0. For the wavelength of A = 40Mm we find Cisobaric ~ — 8 x 
10“® in individual cells for certain snapshots, without clear correspondence 
to the evolution of the blobs. 

Concluding, our results agree with previous studies on the fact that the 
isobaric thermal instability criterion depends strongly on perturbation wave¬ 
length and seems less appropriate to explain catastrophic cooling during 
condensation formation in solar corona. The isochoric instability threshold 
gives a rather nice correspondence with blob formation regions. Hence for 
our case, the isochoric thermal instability criterion is the most appropriate 
to examine blob formation. 

6.3. Hydrodynamic buoyancy and Atwood numbers 

If for whatever reason we arrive at a configuration where heavier fluid is 
positioned above lighter fluid in a gravitational field, this configuration can 
be gravitationally unstable leading to Rayleigh-Taylor (RT) development in 
pure hydro setups. The system will search for a more stable state by forcing 
the heavy fluid to change position with the lighter one and fall down. The 
classical Rayleigh-Taylor instability results when small perturbations are ap¬ 
plied at the interface between the two fluids and prevent the interface from 
keeping a perfectly flat shape. A small perturbation at the interface then 
grows exponentially with a growth rate of: 


exp(A^t), with N = and A = 

Pheavy^ Plight 

where N, is the growth rate, k is the spatial wavenumber and A is the 
Atwood number indicating the density disparity in a gravity field quantified 
by g (Chandrasekhar, 1961; Glimm et ah, 2001). In pure Eulerian fluids, 
small wavelengths grow first. 

This simple criterium is drastically modified in MHD, but it is known that 
the density disparity A between the fluids, which takes values 0 < A < 1, 
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Figure 12: Estimation of the unstable region from an isochoric thermal instability criterion 
shown at three times (top to bottom) and for two visualising angles (left-right column). In 
the panels of the left column we demonstrate the plasma densities with coloured isodensity 
surfaces as indicated by the colourbar. The yellow isosurface visualised here in both left 
(with the density isosurfaces) and right columns (on its own) corresponds to growth rates 
above 0.008Hz and captures well the unstable region, as it follows the changing location 
of the condensed material in the corona with time in our simulation and highlights the 
blobs’ birth places. 
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determines the morphology and the evolution of the plasma undergoing non¬ 
linear Rayleigh-Taylor instability, since for A close to 0, Rayleigh-Taylor 
instability shows a symmetric mixing for both finger-like plasma structures 
for the heavy down falling fluid and bubble-like ones for the light fluid that 
rises; while for A close to 1, falling spikes have larger growth rates and 
penetrate deeper in the opposite region than rising bubbles (Mikaelian, 2014). 



0 0.1998 0.3996 0.5994 0.7992 0.999 

Atwood number 

Figure 13: Atwood number distribution for all the blobs throughout the entire simulation. 
The vast majority of the blobs appears to have Atwood numbers close to unity, suggesting 
a blob density of about 20 times higher than the plasma elements underneath the blob. 


In order to quantify the Atwood number for each individual blob in our 
simulation, we set as Pheavy the maximum density of the blob and search 
underneath the blob on the same column (t/—direction) starting from the 
first grid cell just underneath the blob and moving downwards, for a grid 
cell with density pught smaller than the density of the centroid. There can 
be grid cells underneath the blob with similar density, but hotter material. 
In our simulation the vast majority of the blobs take on Atwood numbers 
of the second limiting case, i.e. values close to unity, as is evident from the 
distribution in figure 13. For the extreme cases we have the following: 

• for A = 0.00145 ^ puob = 1-0029 • Pcorona, 
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• for jA 0.99 Pblob ^ ^9 * Pcorona’ 

As already mentioned, the Atwood number examines Rayleigh-Taylor in¬ 
stability from the purely hydrodynamic perspective, quantifying the density 
disparity between the plasma elements of interest and their surroundings and 
particularly the layer just underneath them, ignoring any influence the mag¬ 
netic field might have on the event. The distribution of the Atwood number 
for the blobs throughout the entire simulation clearly shows that most blobs 
get Atwood numbers very close to unity, which suggests that the conden¬ 
sation centres have about 20 times higher density than their surroundings. 
On the other side of the distribution range we conclude that those blobs 
with small Atwood numbers close to 0 have a density only slightly higher 
than their surroundings. We now address their mixing tendency according 
to more relevant magnetohydrodynamic criteria. 

h./^. Interchange instability and Brunt-Vdisald frequencies 

After the blobs are created, they are situated in low ft regions, so gravita¬ 
tional stratification can cause interchange instability to take place and drive 
their evolution. Thermal instability may still control the delicate growth 
process of individual blobs and sympathetic runaway cooling. Hence, new 
blobs continue appearing at new regions due to catastrophic cooling. After 
their creation, blobs can get born in a gravitationally unstable situation and 
then they start moving to help the configuration find a new more stable state. 
When gravity projected along field lines wins, they fall back to the transition 
region. 

To examine the role played by interchange instability, likely causing the 
blob motion after they are formed, we quantify the relevant Brunt-Vaisala 
frequencies. These are known to appear in standard instability criteria for 
gravity driven modes in up to strong magnetic fields. 

We will use again various instability criteria to trace and follow the evolu¬ 
tion of the blobs and their circulation in the corona. We will find that buoy¬ 
ancy frequencies trace the unstable regions very well and follow the plasma 
condensations as they move around the corona from the moment they are 
formed onwards. We use different quantification strategies for two relevant 
Brunt-Vaisala frequencies. 

The first method in a sense ignores the magnetic topology of the con¬ 
figuration and only takes into consideration the vertical components of the 
gradients (due to gravitational stratification) of density and pressure. This is 
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appropriate in true ID plane-parallel atmospheres, such as an exponentially 
stratified plasma with constant sound and Alfven speed. We then obtain 
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B 


n: 


1 dp f dp 
dy \dy 
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'yp + 3“^ dy ) ’ 


(19) 

( 20 ) 


where is the classic Brunt-Vaisala frequency and is the classic mag¬ 
netically modified Brunt-Vaisala frequency. The latter is especially relevant 
for perpendicular (k ± B) wave modes. 

The second method tries to do justice to the 3D magnetic topology. In 
translationally symmetric (2.5D) Grad-Shafranov type equilibria between 
gravity, Lorentz free and pressure gradient, a stability criterium for con¬ 
vective continuum instability (CCI) modes uses projections on the field lines 
(Blokland & Keppens, 2011). Since we do not have translational symmetry, 
and we want to avoid the details of a straight field line (SFL) representation, 
we arrive at estimations of the Brunt-Vaisala frequencies, taking into account 
the three-dimensional configuration of the magnetic structure of our system 
by 
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( 22 ) 


where ^be Brunt-Vaisala frequencies and ^be magnetically 

modified Brunt-Vaisala frequency, both projected on the magnetic field lines. 
In the actual CCI criteria, one must use the (SFL) poloidal field and project 
along a flux surface. 

We used three different approximations to projected Brunt-Vaisala fre¬ 
quencies, that differ in the way we include the magnetic field influence in our 
calculation. The three different approaches take Bp in the above formulae as 

1. Bp = By, S^) 

2 . Bp = {B^,By,Pi) 

3. Bp = (0, By, 0) 
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Figure 14: Here, we demonstrate for the snapshot that corresponds to 235min of physical 
time the Brunt-Vaisala frequencies, projected on field lines. We show frequencies ^ 
(brown) and (purple) as isosurfaces in the left and right columns, respectively. Ad¬ 

ditionally, in all the panels we visualise in coloured isosurfaces the densities as indicated 
by the colourbar to capture all the regions of high density, i.e. the solar atmospheric 
regions underneath and the blobs inside the corona. We present quantifications done with 
three different ways to approximate field projections, with from top to bottom using a) 
Bp = {Bx, By, Bz), b) Bp = {Bx, By, 0), c) Bp = (0, By, 0). 
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Even slightly negative values of the squared Brunt-Vaisala frequencies, as 
quantihed according to the equations above, indicate unstable states. How¬ 
ever, we explained earlier that due to the resolution of our simulation and 
our time resolution for data saves, we only care for instabilities with growth 
rates corresponding to time intervals of about 1.43min, i.e. the time interval 
between two sequential snapshots. 

The comparison is presented in hgure 14 at time t = 235min. From that 
hgure, we can conclude that the most appropriate method for tracing the 
unstable regions using projected Brunt-Vaisala frequencies is the one that 
corresponds to the total three-dimensional magnetic held, namely the top 
panels in hgure 14. This takes into account the three-dimensional variations 
of density and pressure. The unstable regions estimated in this way surround 
the blobs throughout their evolution (when we make animated views for all 
times) and appear to be a very consistent method to trace the important 
areas of the simulation domain, where the dynamical evolution is critical for 
the condensation centres. The second row, where we use only x, y— variation 
(in line with the initial z—invariance of the setup) of density and pressure, 
overestimates the unstable regions, indicating the whole area underneath the 
loop where the blobs are located. The bottom row is least appropriate, as it 
only takes into account the vertical component of the density and pressure 
gradient and therefore does a bad job localising the perturbed regions. In¬ 
stead, it indicates that the parts of the loops closer to the footpoints are most 
unstable, which is not true, as we observe no blobs or interchange dynamics 
there at all, at that time of the simulation. Hence, the full 3D prescription 
of ah three quantities, i.e. the magnetic held and the 3D gradient of density 
and pressure, is the most accurate method to estimate the unstable regions 
to MHD instability, using both and frequencies. 

The held line projected Brunt-Vaisala frequencies conhrm that we have an 
overall stable conhguration with well dehned unstable regions that surround 
the blobs, throughout their circulation in the corona. They point out as the 
unstable region the top part of the internal loops that form the quadrupolar 
arcade system. We are led to conclude that the blob dynamics is driven by 
some variant of the CCI development. 

On the other hand, the classic way of estimating the Brunt-Vaisala fre¬ 
quencies, taking into account only the vertical variation of pressure and den¬ 
sity, rather than their three-dimensional gradients in the physical domain, 
estimates a larger unstable region, that extends in a large portion of our 
simulation box. This is shown in hgure 15 for a time t = 241min, where we 
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Figure 15: In the four plots above we visualised for the snapshot at about 241min the 
results of four different criteria to estimate the unstable regions in our simulation box as 
they were presented in the text for both Brunt-Vaisala frequencies, N‘^ (left column) and 
N‘^ (right column). We observe that the magnetohydrodynamic approach, i.e. the full- 
held projected Brunt-Vaisala frequencies (bottom row) on the magnetic held lines, nicely 
correlates with the blob regions. The classic approach (top row) that takes into account 
only the vertical component of the gradient of the density and the pressure identihes also 
the (higher P) regions at the domain top, which also seem to show structures as in hgure 
6 . 


33 







contrast it with the full-field projected quantification in the bottom panels. 
More specifically, using the classic formula two unstable regions emerge. One 
on the lower part of the box at the magnetic dip region under the locations of 
blob formation following also the blobs’ evolution, as well as another region 
on top of the simulation box at the upper part of the loops of our magnetic 
topology. They correspond to the dynamics seen at the top of our arcade 
system, also seen earlier in Figure 6, where the plasma beta is higher. 

Still, the most relevant criterion to estimate the unstable regions in low (3 
regions is the one using the three-dimensional gradients of density and pres¬ 
sure and the three-dimensional magnetic field, that captures quite accurately 
the region surrounding the blobs and highlighting the dynamically perturbed 
areas of the domain. 

In figure 16, we demonstrate how the full-field projected Brunt-Vaisala 
frequency-criterion captures the condensation regions throughout the 

simulation following the blobs as they circulate in the solar corona. We do 
so by using six snapshots corresponding to different physical times from the 
moment we observe the first blobs (about 205 min of physical time for the 
top left panel) till the final stages of the simulation (254 min for the bottom 
right panel of figure 16). 

6.5. Additional wave patterns 

There are several different physical phenomena and dynamic processes 
that take place inside our physical domain. As expected and in accordance 
with previous studies of prominence formation (Keppens &: Xia, 2014) we 
see several wave patterns all over the simulation box. The representation of 
the Lorentz force magnitude indicates most evidently all the magnetohydro¬ 
dynamic interactions, as demonstrated in figure 17. However, many of the 
observed wavy features are just an artefact of our (closed) boundary condi¬ 
tions, with repeated reflections at our boundaries as the simulation proceeds. 
Still, we are interested in localised waves that have an actual physical trigger 
at some point in the simulation and such features we will now try to describe. 

From the plethora of wave features that are shown in the visualisation of 
the Lorentz force magnitude all over the simulation domain (figure 17), we 
noticed some clear wave patterns appearing in profiles of other quantities as 
well. One very characteristic example is the region on the top right part of 
the arcade in this plane x = z, as shown with the white line in figure 17. 

There, a moving perturbation seemed to be localised, so we went on to 
quantify the evolution along this line. Our goal is to calculate the propagation 
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Figure 16: We present six snapshots at different times for the unstable region as it is 
estimated by the full-field projected Brunt-Vaisala frequency-criterion A/"^ With brown 
colour and low opacity we can see the isosurface that corresponds to the unstable region 
and with saturated colour bar the density of the blobs. We conclude that there is a 
satisfying correspondence between the criterion and the position of the blobs at each time, 
as the unstable isosurface always surrounds the condensed matter. 
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Figure 17: The Lorentz force magnitude is demonstrated in the slice that corresponds to 
X = z and cuts the domain diagonally at the snapshot that corresponds to about 205min 
of physical time. Several wave patterns appear throughout our simulation box. The white 
line represents the line along which we estimate the physical variations further on. 
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Figure 18: The wave fronts identified by a Lorentz force magnitude view stacked in time 
along the line shown in figure 17 to reveal their propagating speeds. 


speed of this wave and for that we quantify the values of the Lorentz force 
magnitude on the chosen line and stacked them over time, as depicted in the 
figure 18. 

These appear at the middle part of the line length (total length is ~ 
4Lunit in figure 18) and from the slope we calculate the propagation speed 
of these waves. We got an estimated speed for this position-time diagram 
of Lorentz force of about 6km/s. This apparent speed is about an order 
of magnitude smaller than the speed of the blobs while circulating around 
the lower parts of the corona and it is also about 2 orders of magnitude 
smaller than the local Alfven speed. Animated views show clear indications 
of pressure perturbations moving in an oblique direction. 

In the figure 19, we show the magnetic and plasma pressure variations, 
Pfiuid = p{t)-p{U) ^'^<^Pmagnetic = -5^(tj))/2, along the fine presented 

in figure 17, where as we use the moment when we first observe blobs in 
our simulation ti ~ 205min and for t we take t = 241min. Note that in this 
region, (3 is in the range [0.4,1.1], as shown in the left panel of figure 20 (for 
a specific snapshot at 241min of physical time), so the local sound speed is 
in the range [232.9, 256.19]km/s, as shown in the right panel of figure 20. 
From time series of similar figures of these spatial variations in pressure, as 
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Figure 19: Here we demonstrate the magnetic (red curve) and plasma (green curve) pres¬ 
sure deviations as they were quantified along the line shown in figure 17, at physical time 
of about 241min. The magnetic and fluid pressures are out of phase. 





in figure 19, we can estimate the order of magnitude of the wavelength of the 
oscillations as well as their period. We estimated a wavelength A of about 
2.5Mm and a period of about 257s, giving a phase speed of the order of 
Vphase ~ lOkm/s. 

We conclude that the magnetic pressure and the plasma pressure are out 
of phase. This is a typical characteristic of slow magneto-acoustic waves in a 
uniform medium. Even though, our medium is far from uniform, this result 
is a clear indication for the characterisation of the examined waves as slow 
magneto-sonic. 
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Figure 20: Here we demonstrate the plasma /? (left panel) and the characteristic plasma 
speeds locally, i.e. sound speed (red curve in the right panel) and Alfven speed (green 
curve in the right panel), as they were quantified along the line shown in figure 17, at 
physical time of about 241min. 


7. Conclusions 

We use a 3D magnetohydrodynamic simnlation and include thermody¬ 
namic gain-loss terms in a solar quadrnpolar arcade system with magnetic 
dip. Due to evaporation from the chromosphere, eventually localised run¬ 
away condensation events occur and blobs of cool material are formed, as 
a result of thermal instability. We use the thermal instability criteria for 
isochoric as well as isobaric cases to show that the isochoric one is the most 
suitable and gives a satisfying estimation of the nnstable regions where con¬ 
densations create blobs. Indications of interchange instability are evident 
and a continuous coronal rain develops. The interchange instability takes 
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over after the condensations make their appearance, as the configuration is 
gravitationally unstable with dense plasma on top of rarefied plasma. Our 
magnetic field dominates on the lower part of the simulation, as indicated 
by the plasma-/? parameter so that one must exploit field projected stability 
analysis. Blobs are heavily affected by the presence of the magnetic field and 
they seem to accumulate in the middle of our configuration, while splitting 
and merging to finally follow the magnetic field topology and fall back to 
the transition region. The latter occurs along the footpoints of the arcades 
with an accelerated speed. The estimated mass captured by the blobs is in 
agreement with the order of magnitude of observed mass loss rate due to 
coronal rain events, i.e. 5 x lO^g according to Antolin & Rouppe van der 
Voort (2012). We examined commonly used quantifications of Brunt-Vaisala 
frequencies, trying to improve the way the magnetic field topology is taken 
into account. They all help to identify the buoyant unstable regions with 
one of them taking into account the three dimensional magnetic field and 
projecting the pressure and density gradient along the field lines. When only 
taking into account the vertical component of pressure and density gradients 
we also find the (top) unstable low field regions in our domain, where no 
blobs are found but where other subtle velocity-temperature fluctuations do 
occur. We experimented with three different variations of the projection on 
the field fines and concluded that overall, the most trustworthy method for 
tracing the unstable regions is the projection with the full three dimensional 
magnetic field topology. This is also consistent with the CCI criterium for 
more idealised setups. 

A plethora of waves becomes obvious throughout the simulation box when 
visualising the Lorentz force variation. Oscillations along the magnetic field 
lines were also observed in intensity variations in O’Shea et al. (2007) sug¬ 
gesting the occurrence of standing waves and propagating perturbations. We 
found rather localised structures on the profiles of other physical quantities 
as well. We analysed a typical case and estimated its wave propagation speed 
to be about 6km/s and a phase speed of the order of lOkm/s. These speeds 
are of the order of magnitude as the redshifts of downflows observed in 1548A 
and 630 A and the emission stops sharply after the condensed material has 
fallen back to the chromosphere (Muller et ah, 2004b). From the pressure 
variation we conclude that they are magneto-acoustic in nature. 

The instability criteria that we used are only strictly applicable in the 
linear regime/phase of the instabilities. In our paper, we show how they 
can also serve to gain some insight on the physical processes, even when 


40 


non-linear phenomena are the ones dominating. 

Concluding, this simulation is one of the hrst ones that captures the mass- 
cycle in a coronal loop system initiating from footpoint-localized heating, 
causing chromospheric evaporation, triggering thermal instability that will 
result in plasma condensation and downflow, in full 3D setups. This physi¬ 
cal sequence of processes is related to phenomena such as prominences, flares 
and coronal rain, that are overall closely intertwined, as indicated in previous 
studies (Ahn et ah, 2014; Antolin et ah, 2012a; Antolin & Verwichte, 2011; 
Karpen et ah, 2001; Kawaguchi, 1970; Keppens & Xia, 2014; Murawski et 
ah, 2011; Oliver et ah, 2014; Shimojo et ah, 2002). A quadrupolar arcade 
system with a dip, a conflguration usually met in prominences, was used in 
this experiment, as a natural continuation of previous simulations in coronal 
rain and prominence related phenomena (Fang et ah, 2013; Keppens & Xia, 
2014). Catastrophic cooling taking place at a height of about 25 Mm, in 
agreement with Antolin & Verwichte (2011), leads to the formation of cool 
plasma condensations, that emit in EUV, as demonstrated in figures 7, 8. 
Cool and dense plasma blobs are observed in cool chromospheric lines, such 
as Ha and Ca II H, as they descend. A wide range of velocities was revealed 
from a few km to almost 100 km s“^, in agreement with previous studies 
(Antolin & Rouppe van der Voort, 2012; Antolin et ah, 2012b, 2010; Fang 
et ah, 2013; Muller et ah, 2004a; Oliver et ah, 2014; Shimojo et ah, 2002). 
Temperatures with peak at 20,000 K and plasma-/? values of maximum 0.557, 
were found, both within ranges already reported (Ahn et ah, 2014; Shimojo 
et ah, 2002), (Antolin & Verwichte, 2011). Regarding our findings on the 
magneto-acoustic wave, Murawski et al. (2011) also found that magneto¬ 
acoustic waves were linked to cool dense descending blobs, while Oliver et al. 
(2014) supports that down falling blobs emit small amplitude sound waves 
of periods of about 100s. It must be stated that our magnetic field strength 
resembles more that of quiet Sun, rather than an active region case. Mag¬ 
netic fields were estimated to an intensity in the range of 8-22G (Antolin 
et ah, 2012a; Antolin & Verwichte, 2011), while at blob length scales, i.e. 
a few hundreds of kilometres, no perpendicular motion to the local mag¬ 
netic field is expected according to Antolin & Rouppe van der Voort (2012). 
Finally, the cool condensed plasma blobs appear as thin and elongated struc¬ 
tures morphologically at the projected view of figure 11 towards the end of 
the simulation, as observed in several cases (Antolin et ah, 2010; Antolin & 
Verwichte, 2011; Antolin et ah, 2012b). 
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